Concentration and mass dependence of transport coefficients and correlation functions in binary 

mixtures with high mass-asymmetry 
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Correlation functions and transport coefficients of self-diffusion and shear viscosity of a binary Lennard- Jones 
mixture with components differing only in their particle mass are studied up to high values of the mass ratio (i, 
including the limiting case fi, = oo, for different mole fractions x. Within a large range of x and (x the product 
of the diffusion coefficient of the heavy species D% and the total shear viscosity of the mixture r\ m is found to 
remain constant, obeying a generalized Stokes-Einstein relation. At high liquid density, large mass ratios lead 
to a pronounced cage effect that is observable in the mean square displacement, the velocity autocorrelation 
function and the van Hove correlation function. 



I. INTRODUCTION 

The dynamic properties of binary fluid systems where one 
particle species differs from the other only in size, mass or 
both of these parameters have been the subject of a large num- 
ber of studies during the last years H H B H & & B 11 i, 
[HI [IH EH] • The increasing interest is, on the one hand, due 
to the fact that such systems serve as simple models for col- 
loids and micellar solutions, which are of prime importance in 
many scientific areas such as biology or biochemistry, on the 
other hand it is sparked by the rapidly growing capabilities 
of modern computer hardware which allows us to investigate 
parameter ranges and system sizes that were not accessible 
before. 

In general, there are two important limiting cases that can 
be focused on. The first and best-investigated so far is the so- 
called 'tracer' or Brownian limit of one single heavy and/or 
large molecule suspended in a solvent of light particles (in- 
finite dilution). This is especially interesting because for the 
case of a macroscopicalry sized and in comparison to the sol- 
vent infinitely heavy tracer particle, there is a simple relation 
between the tracer diffusivity D and the viscosity of the sol- 
vent tj. This Stokes-Einstein (SE) relation states that 

k B T 



D = 



where ks denotes Boltzmann's constant, T the temperature, 
and C is a numeric constant depending on geometric bound- 
ary conditions. Several studies have been devoted to the 
Brownian limit |l|, |2|, |4[ la, la ll ill , some of them especially 
to the question, for which range of tracer mass and size this 
purely hydrodynamic relation also holds on the microscopic 
level ||5l \M. Depending on whether the mass is changed at 
constant size ratio or not, the SE relation was found to hold 
for mass ratios larger than 10 |Ht] and larger than 100 10], re- 
spectively. Above these values, the tracer diffusion was con- 
sidered mass-independent. (In this context, it is also interest- 
ing to note that even for a pure simple fluid the SE relation 
was found to hold in a large part of the phase diagram JUt].) 

Considerably less studies exist on the approach to the 
Brownian limit (small but finite concentration of heavy par- 
ticles, JH HI [l0[]) and mixtures with larger mole fractions of 
the heavy component (IT 



The second limiting case corresponds to the scenario when 
the mass of the heavy species goes to infinity. For the case 
of a single Brownian particle, this situation was covered in 
y|, where it is explained why it makes sense to attribute a 
finite friction coefficient to an immobile particle. For a fi- 
nite concentration of heavy particles, on the other hand, the 
infinite-mass limit effectively transforms the system into a 
one-component fluid in a random porous matrix of fixed ob- 
stacles that takes up a finite fraction of volume. As a conse- 
quence, there exists a percolation transition at large density 
and concentration of the heavy component which ultimately 
prevents the light particles from diffusing through the system. 

But also in the case of finite mass ratio and concentration, 
the heavy particles influence the dynamics of the light ones. 
Since different masses induce different time scales, the heavy 
molecules act as a cage for the light ones, stalling their dif- 
fusion until they themselves have finally moved significantly 
from their starting position. This leads to so-called hopping 
processes, characterized by particles moving from one cage 
to the next. These complicated dynamic processes, which also 
occur near the glass transition in supercooled liquids lfl4l[l5ll . 
have made it difficult to tackle such systems by theory, which 
is why many of the existing studies rely mostly on computer 
simulations. The most successful theory at present is the 
mode-coupling theory (MCT), which in its general form in- 
corporates a mathematical description of hopping processes, 
but also in its idealized form yields already rather good results 
for glassy systems ifTilflr^ fr AfTsIl as well as highly mass- or 
size-asymmetric mixtures [12]. Recently, the mean-field the- 
ory of Tokuyama has also been applied successfully to diverse 
glass-forming systems lfl9l I20I I21I l22l f23n . 

Now, the goal of our work is to study the approach to the 
infinite-mass limit in an asymmetric binary mixture at finite 
concentrations of the heavy component, using molecular dy- 
namics calculations on as simple a model as possible. In par- 
ticular, we have chosen a truncated and shifted Lennard-Jones 
interaction potential with a cut-off radius of r c = 2.5, where 
both species have equal interaction strength and particle size. 
We perform simulations at two state points away from the crit- 
ical region, one with moderate and one with high liquid den- 
sity (the phase diagram for this system can be found in [ 24]). 
The only two remaining system parameters, the mass ratio /j, 
and mole fraction x, do not influence the static properties of 
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the system at all. What we are interested in are the transport 
coefficients of self-diffusion Di and D 2 and the shear vis- 
cosity 7] m of the mixture, especially the dynamics for high /i 
including the case fj, = oo. In order to avoid dealing with 
the percolation th reshold, we focus on small concentrations 
x < 0.2. Although there are no critical fluctuations, due to 
the periodicity of the simulation box the diffusivity is afflicted 
with a finite-size effect that has to be accounted for l25ll . We 
also study a possible relation between D 2 and r\ m (general- 
ized SE relation) and the influence of high /i on the dynamics 
of the light particles (cage effect). 

The paper is organized as follows: In the next section, we 
discuss in more detail the model system and the simulation 
methods we apply. After that, we give a brief overview on the 
basic formulas for the dynamic quantities we are going to look 
at, and subsequently we present the results we have obtained. 
Finally, we end with a conclusion and a short outlook. 



n. MODEL AND SIMULATION DETAILS 

We consider a two-component mixture consisting of N = 
Ni +N 2 particles in a cubic volume V with periodic boundary 
conditions. A pair of particles i, j of species a and j3, respec- 
tively, separated by a distance r = |r$ — | , interacts via a 
truncated and shifted Lennard- Jones potential <f> (r), given by 



^ ( r ) = j ( t >LJ M ~ 4> LJ ( rc ) ' r <r c 



r > r c 



and 



4>lj (r) = 4£ Q/3 



(2) 



where r c = 2.5, and the two species have the same inter- 
action strengths e\\ = e 22 = £12 = s and particle sizes 
an = (T22 = C12 = c, but different masses ra 2 > mi. From 
the point of view of statics, such a system is identical to a one- 
component Lennard- Jones fluid, the dynamic properties, how- 
ever, will of course depend on the mass ratio /i = 777,2 /mi and 
the concentration specified by the mole fraction x = N2/N. 
In the extreme limit of /1 — > 00 we are effectively dealing with 
a system of Ni particles moving in a disordered matrix of A<2 
fixed obstacles of the same size. 

We performed MD simulations in the TWT-ensemble us- 
ing a Nose-Hoover thermostat [26, |27l H^] in the formu- 
lation of Martyna et al. 13011 and a velocityVerlet integration 
scheme. A multiple time step algorithm 113 ill was applied in 
order to deal with the different time scales due to the high 
mass asymmetry of the two mixture components. Simulations 
usually lasted at least 2 x 10 6 time steps, where one time step 
of the light species was chosen as 5t — 0.005y/mia 2 /e, with 
equilibration times of 2 x 10 5 time steps. 

Simulations for infinite mass ratio /1 = 00 were performed 
by fixing the positions of the heavy particles at random con- 
figurations obtained from short simulation runs under identi- 
cal thermodynamic conditions but with equal masses for both 
species. Usually, the reported simulation results were aver- 
aged over 10 different configurations, which turned out to be 
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FIG. 1: (Color online) Points: MD data of the MSD of the light 
species divided by time for x — 0.2 and mass ratios /1 = 2 (red), 
(i — 20 (purple), n — 500 (blue) and fi = 00 (green). Solid lines: 
Mean field results, fitted for the mean free path length which is 
0.15, 0.15, 0.155 and 0.164, respectively. Dashed lines: diffusion 
coefficients obtained from the Green-Kubo relation, multiplied by 6. 
The density is p = 0.9 and temperature T = 1. 



a large enough number, since we found the dependence on the 
exact configurations to be quite small. Also in the case of fi- 
nite mass ratios each result is obtained as an average of several 
simulation runs in order to improve the statistics and calculate 
standard deviations of the dynamic quantities. 



The question how the limit /1 — > 00 is to be interpreted 
is a delicate one [1]. Since the heavy particles have infinite 
mass and zero velocity, their momentum P2 = m 2^i 
is in principle undefined. However, the total momentum of 
the light particles Pi = 53i=i m i v « i s known at any time, 
and therefore by requiring the total momentum P of the sys- 
tem consisting of light and heavy particles to be fixed and (by 
convention) equal to zero - just the same as in the simula- 
tions with finite /i - we can assign them the finite momentum 
P2 (t) = —Pi (t) . In this way, we define the system at any 
time step t as the limiting case of a system where 7772 goes 
to infinity and v; goes to zero for all i S {1, . . . , A^}, while 
P2 is held constant at the value — Pi (7;) . The same procedure 
was applied in (jl to calculate the momentum autocorrelation 
function of a single Brownian particle with infinite mass. 



In the following, all quantities will be given in reduced 
units: energy is measured in e, length in a, mass in 777,1 and 
time in -J niia 2 je. Temperatures are given in e/kg, densities 
inl/cr 3 . 
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III. THEORETICAL BACKGROUND 

A. Mean square displacement and diffusion coefficients 

The mean square displacement (MSD) of a particle of 
species a is defined as 



-, N a 

(Arl (t))= w Y,([r i {t)- ri (0)] 

l — l 



(3) 



where (t) is the three-dimensional trajectory of particle i, 
and (•) denotes a canonical average. In the case of normal 
diffusion according to Fick's law, the MSD obeys the so- 
called Einstein-Helfand relation ll32ll . which states that the 
self-diffusion coefficient of light or heavy particles is given 
by its slope at large times according to 



D a = lim 

t — >oo 



<Ar^ (<)) 



(if 



(4) 



Another way to obtain the diffusion coefficients of the two 
components is via their velocity autocorrelation functions 
(VACFs) ip a (t), defined as 



4>o 



1 N a 

«(*) = 3]^E< v *(*) v < (°)) 



(5) 



i=i 



These allow us to calculate the self-diffusion through the 
Green-Kubo relation 1331 



1p a (t) dt. 



(6) 
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FIG. 2: (Color online) Normalized velocity auto-correlation func- 
tion il>i(t) of the light (a) and ^(t) of the heavy particles (b) in a 
system with T = 1, p = 0.9, s = 0.2 and different mass ratios ^i. 



B. Mean Field Theory 



C. Van Hove correlation function 



Tokuyama has established a mean field theory 11911 for the 
MSD near the glass transition in colloidal suspensions and 
molecular systems. In three dimensions, the MSD denoted 
by M (t) is described by the nonlinear differential equation 



j t M (t) = 6D + 6 [v 2 t - D] e - M W 
with the formal solution 

M (t) = QDt + l 2 \n\e- 6Dt/l2 



(1) 



6D 2 



1 - (l + QDt/l 2 



-6Dt/l 2 



(8) 



Here, D is the self-diffusion coefficient, I is the mean free 
path, and vo denotes the average velocity of an atom. We 
have applied this theory to the MSD data for the light particles 
from our MD simulations, taking the value of the diffusion 
coefficient from the Green-Kubo results, and fitting the value 
of / in Eq. (© to the MD data. 



Space-time correlations between two particles in a pure 
fluid system are described by the van Hove correlation func- 
tion (VHCF) G (r, t) , which is defined as Q 



G(r,i) 



+ r, (0) - Tj (t)} 



(9) 



For an isotropic system, G (r, t) d 3 r = 4irr 2 G (r, t) dr gives 
the probability to find a particle at time t a distance r from the 
origin, provided that at time t = a particle was located at the 
origin. The VHCF can be separated into a self- and a distinct 
part, 



G(r,t)=G s {r,t)+G d (r,t), 



(10) 



where G s includes the terms with i = j, and Gd those with 
i j. The self-part, on the one hand, is the time-dependent 
conditional probability density that a particle moves a distance 
r = |r (0) — r (t)| during time t.Att = 0, G s (r,0) = S(r), 
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FIG. 3: (Color online) Non-Gaussian parameter 02(f) of the light 
species for a concentration of x — 0.2 and different mass ratios. The 
density is p = 0.9 and temperature T = 1. 




(a) (b) 

FIG. 5: (Color online) Normalized distinct part of the VHCF, 
G]i(r,t)/p, as a function of r at constant time t, with x — 0.2 
and ft = 10 4 (a) and p = 00 (b). The times are t = T x 0.015, with 
i given by the number next to each curve. The density is p — 0.9 and 
temperature T = 1 in both cases. 



diffusion coefficient of the system, 




G s (r,f) 



(4ttL><) 



.3/2 



exp 



r 2 \ 
4~Dt J 



(12) 



The distinct part, on the other hand, represents the con- 
ditional probability density of finding a particle at time 
t a distance r apart from the location of another parti- 
cle at time t = 0. At t = 0, Ga (r, 0) = pg (r), and 
lim r _ >oo Gd (r, t) = lim^oo Gd (r, t) = p. The normaliza- 
tion is 4-7T / r 2 G d (r, t) dr = N — 1. 

In a binary mixture one has to differentiate between the dif- 
ferent species a, and the corresponding VHCFs are defined 
as 

G a s {T,t) = ^(Y d 5[v + T i {0)-v i (t)]\, 



and 



GT (r,t) 



iV Q (JV Q - 1) 



N a 



^<5[r + r, (O)-rj(t)] 



(13) 



FIG. 4: (Color online) Normalized distinct part of the van Hove 
function G^ 1 (r, t) / p. The concentration of the heavy particles is 
x = 0.2, and the mass ratio is p — 10 4 (a), p = 3 x 10 s (b) and 
p = 00 (c). The density is p = 0.9 and temperature T = 1 in all 
cases. 



^(M) = ^(ED^M0)- rj (t)l 



(14) 



whereas Hindoo G s (r,t) = lim i _ ) . o0 G s (r,t) = 1/V « 0. 
It is normalized to unity, 4-7T J r 2 G s (r,t)dr = 1, and con- 
nected to the MSD via the relation 



(Ar 2 (t)) = / r 2 G s (r,t)d 3 r. 



(11) 



For large enough r and t, the self -part approaches a Gaussian 
distribution whose width grows with \Dt, where D is the 



D. Shear viscosity 

The total stress tensor of a one- or multicomponent system 
is given by ll35ll 



N 
i=l 



N r x r V 

ij ij 



E 

j>l 



0' (Tij) 



(15) 
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Now, the stress tensor of component a in a mixture can be 
defined as 



E 



m i v i v i 



N x r V 
ij ij 

j>i 



(16) 



such that the total stress tensor is the sum of those of the two 
species, 



a 



a 



■ xy 



(17) 



Then one can write the stress-stress auto- and cross- 
correlation functions j] a p (t) with a, [3 E {1, 2} as 

»?«/»(*) = <<C(*)o£y(0)>, ( 18 > 

and thus the total correlation function is given by 

rj (t) = (<j xy (t) a xy (0)) = mi (t) + V22 (t) + 2 m2 (f) . (19) 

The corresponding shear viscosities are computed via the 
Green-Kubo formula 



1 



Vaf3 



Vk B T Jo 



r\ a p (t) dt, 



(20) 



and the total shear viscosity of the mixture rj m can be obtained 

as 



J/m = m + m + 2r? 



12; 



(21) 



where we write m instead of mi an( l V2 instead of r)22- The 
reason for separating r/ m into two contributions is that by do- 
ing so we can compare the results for finite /i with those for 
fi — oo. In this limit, the total viscosity goes to infinity while 
?/i stays finite and approaches m (/•* = oo), the viscosity of the 
light particles moving through the porous matrix. 

The numerical procedure we used to calculate Eq. (|20]) 
was to store the values of a" y (t) , a yz (t) and <r" z (t) at ev- 
ery third time step on disk, and perform the integration via 
a Fast Fourier Transformation after the simulation, averaging 
over the three tensor elements in order to decrease statistical 
errors. In this way, it was not necessary to know the maximum 
integration time beforehand. 



E. Stokes-Einstein relation 

If we consider a macroscopic sphere ('tracer') of radius R 
moving at a velocity V in a liquid with shear viscosity rj, ac- 
cording to Stokes' law, the frictional force acting on the sphere 
is given by 



F = -CV, 
where £ denotes the friction coefficient, 

C = Cnr/R, 



(22) 



(23) 



and the constant C depends on the boundary conditions. In 
the case that the viscous fluid sticks perfectly to the surface 
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FIG. 6: (Color online) Normalized distinct part of the van Hove 
function G^ 2 (r, t) / p. The concentration of the heavy particles is 
x — 0.2, and the mass ratio is /i = 10 4 (a), fi = 3 x 10 s (b) and 
jU = oo (c). The density is p = 0.9 and temperature T = 1 in all 
cases. 



of the sphere (rough surface; stick boundary condition), i. e. 
the fluid velocity is v = V everywhere on the surface, C is 
equal to 6. If, on the other hand, one assumes that the fluid 
slips perfectly over the sphere (smooth surface; slip boundary 
condition), the value obtained for C is 4. In this case, only the 
normal component of the fluid velocity at the surface of the 
sphere is equal to that of the sphere velocity (v± — V±; no 
fluid can enter or leave the sphere), and the tangential force 
at the surface is zero. These two values of C can be found 
through purely hydrodynamic calculations, see for example 



Now, according to Einstein ll38ll . the friction coefficient £ 
is inversely proportional to the diffusion coefficient D of the 
sphere with the thermal energy fc^T as the constant of pro- 
portionality, 



D 



k R T 



(24) 



Equation ( 124b is known as the Stokes-Einstein (SE) relation. 
In the Brownian limit (tracer limit) of a binary isotopic mix- 
ture (only one particle of the heavy component), the Stokes- 
Einstein relation was also verified to hold in the form J1101 



D B = 



k R T 



CittisRh 



(25) 



where Dg is the diffusion coefficient of the Brownian particle, 
tjs is the shear viscosity of the solvent (light component), and 
Rh is the so-called hydrodynamic radius. 
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FIG. 7: (Color online) Self-part of the van Hove function for the 
light particles Gl(r,t) times 4ivr 2 . The concentration of the heavy 
particles is x = 0.2, and the mass ratio is p = 10 4 (a), p = 3 x 10° 
(b) and p = oo [(c), (d)]. The density is p = 0.9 and temperature 
T = 1 in all cases. Low values are shaded in blue, high values in 
red. The plot range in ^-direction is [0,2], the difference between 
two contour lines is 0.12. Clipped regions are shown in white. The 
thick black lines correspond to the root mean square displacement 
\J (Arf (t)}. The curves shown in (d) correspond to the times t = 
2 l x 0.015, with i given by the number next to each curve. 



IV. RESULTS 

A. Mean square displacement and cage effect 

We first want to focus our attention to the diffusion of the 
light particles. Figure Q] shows the MSD of species 1 divided 
by t in mixtures with fi = 2, fj, = 500 and /i = oo (circles) 
for a density of p = 0.9, concentration x = 0.2 and tem- 
perature T = 1, along with the fitted MF curves [Eq. ([8]); 
solid lines]. The dashed horizontal lines indicate the values 
of the diffusion coefficients times 6, obtained from indepen- 
dent calculations via the Green-Kubo formula ©. Obviously, 
there is consistency between the two calculation routes in the 
cases of finite fi, whereas for infinite mass ratio the MSD does 
not reach a linear-time behavior within the length of the sim- 
ulation. In general, we can distinguish three regimes of the 
MSD: the quadratic regime for small times where the parti- 
cles move ballistically at constant velocity, the linear regime 
with the usual diffusion as given in Eq. © for large times, and 
an intermediate region of anomalous diffusion. Such a diffu- 




FIG. 8: (a) Trajectory of a particle of the light species in a system 
with fixed heavy particles (p = oo), mole fraction x = 0.2, density 
p — 0.9 and temperature T = 1. (b) displacement Ar(t) of the same 
particle. 



sion is often explained by the so-called 'cage effect' denoting 
the fact that particles are trapped inside a cage formed by their 
surrounding neighbors for some time, before they can escape 
and diffuse in the usual way. As seen from Fig. 1 the region of 
anomalous diffusion increases with /i. This supports the idea 
13911 that the trajectories of the light particles for large enough 
fi change from relatively smooth ones (Gaussian-like process) 
to intermittent ones with a large amplitude of displacements 
(highly non-Fickian process or activated hopping). This cage 
effect can be seen in several other quantities apart from the 
MSD. 

For example, the VACF of the light particles shows a dis- 
tinct negative minimum, that becomes more pronounced if fi 
increases, as can be seen in Fig. [2] (a). In the VACF of the 
heavy particles, on the other hand, the minimum vanishes with 
increasing mass ratio [Fig. 12(b)] . The position of the mini- 
mum i m j n offers a way to estimate the typical size of a cage. 
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FIG. 9: Part of the trajectory shown in Fig.[8]for t between 1200 and 
1350. The particle is trapped in a potential energy minimum. The 
isosurface corresponds to the total potential of all heavy particles. 
Inside the surface, the potential is lower than -5. 



With a mean thermal velocity of t>o = V3 at T = 1 the light 
particles on average travel a distance vot m \ n = 0.26 until they 
are reflected by the surrounding cage. Together with the par- 
ticle radius 0.5 this results in a cage diameter of about 1.5. 

Another indicator for the cage effect is the so-called non- 
Gaussian parameter (NGP) a 2 , considered as a measure of 
dynamic heterogeneity on intermediate time scales. Based on 
Eq. (fTTT i and the fact that the self-part of the VHCF takes a 
Gaussian form in the case of normal diffusion [Eq. dTzb l. one 
defines [40] 



3(r 4 (t)) 
a 2 (t)= K V >' - 1. 
5 (r 2 (i)) 2 



(26) 



Typically, a 2 takes a value of ~ 0.1 — 0.2 in the normal fluid 
regime 14111 . Our simulation results of the NGP for the same 
system as in Fig. Q]are shown in Fig. [3] As one can see, with 
increasing mass the peak value of the NGP reaches almost 
0.6, starting from the value 0. 12 for fi = 1, indicating that the 
cage effect becomes more pronounced for larger p. This is in 
accordance with the behavior of the MSD seen in Fig. Q] It is 
also visible that for infinite mass ot 2 if) is still nonzero at the 
largest times considered in our simulations. 

Finally, the cage effect should be observable via the VHCF. 
For normal diffusion, the value of the distinct VHCF at the 
origin goes monotonically from to p with increasing time 
IB4I1 . An additional peak appearing at r = indicates that a 
particle has 'hopped' into the cage where the reference parti- 
cle has been at t = 0, which by this time has escaped this cage 
formed by the surrounding particles [14]. Figures [4] - [6] show 
the normalized distinct VHCFs G] l 1 (r,t)/p and G 1 d 2 (r,t)/p 
for the same density, concentration and temperature as in Figs. 
[T]|3]and large mass ratios p = 10 4 ,3 x 10 5 and oo, where 
this effect is well seen. At p = 10 4 [Figs. @](a) and|5](a)], 



G 1 ^ exhibits a peak at r = with about the size of the first 
peak of the pair distribution function g(r). With increasing 
mass ratio, both the height and the width (in time-direction) 
of the peak grows, until for infinite mass ratio [Figs. H](c) an d 
(b)] the peak persists even up to the largest time, which is 
0.015 x 2 19 . In G d 2 , shown in Fig. [6] the effect of increas- 
ing mass is apparent in the fact that the structure of the pair 
distribution function g(r) which is present at t = dissolves 
at later and later times, until it remains unchanged throughout 
the whole simulation time for p = oo [Fig. |6](c)], reflecting 
the 'frozen' configuration of the heavy particles. 

But also the self-part of the van Hove function can reveal 
something about the cage effect and dynamic heterogeneity. 
As we already mentioned, the usual shape for G s (r, t) is a 
Gaussian distribution in r for any large enough time, with its 
width increasing like \f~Dt. This reflects a Fickian process. 
However, if the dynamics is strongly intermittent, it becomes 
heavily non-Gaussian at intermediate time range. For such an 
anomalous or hopping diffusion, the appearance of an addi- 
tional peak is typical [ 14] . We show the self part of the VHCF 
of the light particles G\ (r, t), multiplied by 47rr 2 to obtain the 
probability density, in Fig. [7] for the same systems as in Fig. 
@] and [6] The formation of a multi-peaked structure when p 
takes large values ranging from 10 4 to infinity is clearly visi- 
ble in Fig. [7] This Figure also demonstrates the three-peaked 
nature of Gl for p = oo. The distance between the peaks in- 
dicates that the typical distance of two cages is about the size 
of a particle. 

In order to examine the hopping behavior and cage entrap- 
ment more closely, we have taken a look at the trajectory of 
a single light particle in the course of a simulation run with 
infinite mass ratio. The mole fraction of the heavy species 
was again x — 0.2, the density p = 0.9, the temperature 
T = 1, and the number of particles used in the simulation was 
N = 500. Figure[8](a) shows the obtained three-dimensional 
path through the simulation box (the fixed heavy particles are 
depicted by gray spheres; periodic boundary conditions ap- 
ply), whereas Fig. [8](b) shows the distance Ar(t) from the 
starting position at t = covered by the tagged particle. Both 
figures demonstrate that the particle is repeatedly trapped at 
some place, oscillating around a position with an amplitude 
smaller than the particle size, before it hops again to some 
other trap. Hence this trajectory demonstrates well the inter- 
mittent process, discussed above. In Fig. [8](a) these traps ap- 
pear as black regions where the trajectory passes many times. 
In order to verify that the traps are indeed minima of the poten- 
tial energy landscape created by the fixed particles, we have 
plotted in Fig. [9] a magnified portion of the trajectory, corre- 
sponding to the time period 1200 < t < 1350, during which 
the particle is trapped according to Fig. [8](b). It is obvious that 
for fi = 1 Ar(i) in Fig. 8 (b) would show the known behav- 
ior of a pure fluid where plateau ranges are practically absent. 
Only for large enough values of the mass ratio and interme- 
diate concentrations time intervals of constant Ar(i) become 
visible. Also shown in the figure is a surface of constant po- 
tential energy (the other mobile particles are not included in 
the calculation). It is apparent that the surface forms a kind of 
bag, and the trajectory lies almost completely on the inside of 
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FIG. 10: (Color online) Diffusion coefficient of the heavy (full sym- 
bols) and light species (open symbols) as a function of the mass ratio 
p for concentrations x = 0.05 (green) x = 0.1 (blue) and x = 0.2 
(red). The dashed curves correspond to the linear model (see text). 
The density is p — 0.6 and temperature T = 1.05. 
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TABLE I: MD results for the kinetic energies 
light and heavy subsystems, for various concentrations, mass ratios 
and system sizes. The temperature is T — 1.05 and the density 
p — 0.6 in all cases. 
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FIG. 11: (Color online) System size dependence of the diffusion 
coefficient of the heavy species for p = 10 and different concentra- 
tions. Solid curves are extrapolations to 1/L — 0, dashed curves 
are only guides to the eye. The density is p = 0.6 and temperature 
T = 1.05. 



it, where the potential energy is smaller than on the outside. 
The dimensions of the portion are roughly a x a x ~, and 
there is no heavy particle inside. 



B. Diffusion coefficients 

Simulation results for the diffusion coefficients of light and 
heavy particles, D\ and D 2 as functions of the mass ratio 
H are presented in Fig. [I0]for three different concentrations 
x = 0.05, 0.1 and 0.2. The temperature is T = 1.05 and 
the density p = 0.6. The results were corrected for finite-size 
effects according to the equation | |25[ |42 



D(L) = D(oo) 



a 

T 



(27) 



where a is a fitting parameter. Figure QT| shows the depen- 
dence of D 2 on the system size L for some concentrations 
x = 0.2, 0.1 and 0.03. In each case the mass ratio is [i = 10. 
As one can see, the curves exhibit an increasing curvature with 
decreasing x, reflecting a departure from the 1 /X-scaling be- 
havior predicted by Eq. d27| >. Extrapolations to 1/L = were 
performed using the two data points with the highest particle 
number in each case. 

Another difficulty that is caused by a small number of parti- 
cles with high mass is that the mean kinetic energies, and thus 
the temperatures, of the light and heavy subsystems may de- 
viate from each other. For a single Brownian particle of mass 
tub, this problem was studied in detail by Nuevo et al. JUl. 
For a total number of particles and a mass m of the solvent 
particles, the mean square momentum of the Brownian parti- 
cle (p 2 B ) will differ from its value in the thermodynamic limit, 
3Titlb, by a factor of 



N - 1 



N — 1 + mB/rn 



(28) 



Equation (|28|) , derived in [1], can be generalized to an arbi- 
trary number A^ of heavy solute particles and N% = N — A^ 
light solvent particles, yielding 



= N 1 + (N 2 -l)[i = l-x+(x-l/N)n { 



N 1 + N 2 n 



1 — X + XfX 



Table U gives some examples of measured kinetic energies of 
the two species, E\™ and E^, compared to the value pre- 
dicted by Eq. fl29l) . The thermodynamic limit value for this 
temperature is E^ in = |T = 1.575. In some cases with low 
concentration and small system size the deviations are found 
to reach up to 10%. 
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FIG. 12: (Color online) Total shear viscosity r\ m as a function of the 
mass ratio p for systems with p — 0.6, T = 1.05, and different mole 
fractions x — 0.2, 0.1 and 0.05. The dotted curves correspond to the 
linear model (see text). 



FIG. 13: (Color online) Total shear viscosity rj m and the contribu- 
tions r/i and 772 of the two components as functions of the mass ratio 
p for a system with p = 0.6, T = 1.05, and mole fractions x — 0.2. 
The black dotted curve correspond to the linear model (see text), the 
blue dashed line indicates the value of 771 obtained for p — 00. 



C. Shear viscosity 

For the same systems as in Fig. [I0]we have calculated the 
shear viscosities of the mixtures via Eqs. ( f20b and (|2TT >. The 
results are presented in Figs. [12] and [13] We compare them 
with a simple linear model assuming the mixture is ideal, and 
therefore the total viscosity 77*^ is given by 

r,%=(l-x)r 1 1 +xr 1 2 , (30) 

where r)i and 77° denote the shear viscosities of the two com- 
ponents in their pure form. Since the viscosity of a pure 
fluid scales with the square root of the mass of its particles, 
7/2 = J~pri\, we have 

% = l+x(^p-l). (3D 
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FIG. 14: (Color online) Normalized stress-stress autocorrelation 
function 772 (t) for the heavy component of the system with x — 0.2 
of Fig. E] 



The dotted lines in Figs. [T2l and [T3l were obtained from Eq. 
[3T1 Agreement with the MD data is in general quite good, 
only for x — 0.2 the model overestimates the real values by 
up to 20%. Figure [13] shows additionally the contributions 771 
and 7/2 of the two mixture components as defined in section 
IIIIDI It can be seen that while rj m and 772 are both increasing 
as y/JI for large mass ratios, 771 is growing only slowly and 
reaches the value obtained for p — 00 (blue dashed line) at 
p = 10 4 . 

In Fig. Q3] we plot several stress-stress autocorrelation 
functions 772 (t) = 7722 (t) [see Eq. ([TST ll for the systems with 
x = 0.2 and various values of the mass ratio p. It is obvious 
that with increasing p the relaxation times grow in a similar 
manner as we observed for the VACF ip2(t)- Consequently, 
the numerical integration of Eq. (|20| i has to be extended up to 
very large times t m . dx ~ 100 in order to reach the plateau value 
of 772. 



D. Stokes-Einstein relation 

In order to look for a relation between the diffusion coef- 
ficient of the heavy particles and the shear viscosity, we plot 
D2 from Fig. [l0]and r\ m from Fig. [T2]in a double-logarithmic 
scale in Fig. [15] (full symbols). We observe that these data 
points lie close to a straight line, which might be represented 
by the equation 

D 2 = A Vm a , (32) 

where A and a are fitting parameters (a similar relation was 
also suggested in lfl2tl ). Indeed, a linear fit leads to the val- 
ues a — 1.07 and A — 0.131. Thus, we propose a Stokes- 
Einstein-like relation with a = 1 which yields a value of 
A = 0.122 (solid line in Fig. [TBI . From such a relation one 
may extract an effective hydrodynamic radius Rh by iden- 
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TABLE II: MD results for the diffusion coefficients and viscosities, 
obtained via the Green-Kubo formulas (f6j and l !20t . for various con- 
centrations and mass ratios. The temperature is T — 1.05 and the 
density p — 0.6 in all cases. 



tifying A = ksT /CttRh according to Eq. (|25]) . Assuming 
slip boundary conditions, we obtain Rh = 0.68, which seems 
to be reasonable for our interaction potential. Similar values 
have also been found before, e. g. in Q9[] . 

Equation (f32l > also allows us to apply the linear ideal- 
mixture model OTb for rj m to the diffusivity D2 . Combining 
the two equations yields 



D 2 (x, n) 



D 



1 



M-l 



(33) 



with B = A/i~ii = 0.156. The curves resulting from Eq. (133]) 
for x = 0.05, 0.1 and 0.2 are shown by the dashed lines in 
Fig. [TO] 



E. Concentration dependence 




FIG. 15: (Color online) Fit of the MD data to Eq. <[32j for con- 
centrations x — 0.2, 0.1 and 0.05 and different mass ratios p = 
1, 10, 100, 500, 3000, 10 4 . The straight line corresponds to an ex- 
ponent of -1, and the fit yields a hydrodynamic radius Rh = 0.68 
(assuming slip boundary conditions). The density is p — 0.6 and 
temperature T = 1.05. 



predicted by the SE-relation (|32|) is too low, which is also 
apparent from Fig. Q3] It seems that for not too small values of 
x, the linear model describes the behavior quite well, while for 
small concentrations the deviations are getting larger. Also, 
from the MD data it is not clear whether D2 approaches the 
same value for any fiasx goes to zero. Since finite-size effects 
increase when approaching the Brownian limit, we could not 
answer this question. 

A clearer picture can be given regarding the concentration 
dependence of the shear viscosity, shown in Fig. [16] (b). For 
x — ► 0, Tj m of course approaches the pure-fluid value rf( of 
the light component, which is 0.786(12) for the chosen den- 
sity and temperature. At small concentrations, the function 
77m (x, p) can be approximated by a linear ansatz, 



1 + kr, (//) X + 



(34) 



with a /i-dependent coefficient k v . Comparison with Eq. (l3~n > 
yields k n (fi) = s/fj, — 1 for the simple linear model. The 
slopes obtained from the simulation data (k v = 19 for // = 
500, k v = 70 for pi = 3000 and k v = 186 for /j = 10000) are 
in qualitative agreement with this assumption. In any case, the 
observed numbers are much larger than the well-known value 
of 2.5 proposed by Einstein for a suspension of solid particles 
in a liquid at small concentrations [38, 4~4ll . 



Finally, we investigated the dependence of D2 and rj m on 
the concentration. Figure [16] (a) shows the diffusivity of the 
heavy component as a function of x for p — 0.6, T = 1.05, 
and pi — 10, 100, 500 and 3000. For comparison, we also 
include the curves predicted by the linear model, namely Eq. 
((331) at fixed values of p. In this case, however, we set B = 
DzipL = 1) = 0.187, since otherwise the value of Di{x — > 0) 



V. CONCLUSION 

We have performed extensive MD simulations of binary 
Lennard-Jones fluids whose components are identical except 
for their mass, such that only dynamic properties like trans- 
port coefficients and time correlation functions change with 
varying mass ratio p, and concentration x of the two species. 
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FIG. 16: (Color online) Dependence of the diffusion coefficient of 
the heavy species (a) and the shear viscosity of the mixture (b) on 
the concentration for different mass ratios. The dotted curves in (a) 
correspond to the linear model (see text), the dashed curves in (b) are 
linear fits. The density is p — 0.6 and temperature T = 1.05. 



We found that especially at high liquid densities and high 
mass ratios the large difference in relaxation times of light 
and heavy particles leads to a pronounced cage effect for the 
light component. It can be observed as an intermediate re- 
gion of anomalous diffusion in the mean-square displacement, 
a large maximum of the non-Gaussian parameter, and addi- 
tional peaks in both the self- and distinct part of the van Hove 
correlation function. When tracing the trajectory of a single 
light particle, it turns out that its motion is characterized by 
hopping between separate local minima of the potential en- 
ergy landscape. Thus our study gives in fact an example that 
a rather stable 'solvent cage' can be formed in mixtures just 
because of a strong mass asymmetry effect. 

Furthermore, we established a generalized Stokes-Einstein 
relation between the diffusion coefficient of the heavy compo- 
nent and the total shear viscosity of the mixture that is valid 
in the whole range of mass ratios and concentrations. In order 
to obtain accurate results, it was necessary to correct for the 
system size dependence of the diffusivity, and to ensure that 
the Green-Kubo integral for the shear viscosity has reached its 
plateau value. 

Mass dependence of both viscosity and diffusivity are ap- 
proximately predicted by a simple linear model assuming an 
ideal mixture behavior. For small concentrations, the shear 
viscosity follows a linear dependence on x with a slope going 
roughly as JJi, whereas for the diffusion coefficient of the 
heavy species due to computational limitations no conclusive 
result could be obtained. 

There are several possible ways of extending the results pre- 
sented here. Apart from studying other transport coefficients 
such as thermal conductivity or mutual diffusion, it would be 
interesting to perform a similar investigation close to the crit- 
ical point of the phase diagram. Furthermore it is planned to 
study the percolation threshold at high concentrations of the 
heavy component, and to leave the hydrodynamic regime and 
look at the wave-vector dependence of the dynamic quantities 
(see, e.g., ll45lo . For the latter problem, calculations have al- 
ready been performed the results of which will be published 
elsewhere. 



In particular, we have studied diffusion coefficient, shear vis- 
cosity, velocity and stress-stress autocorrelation functions, 
the van Hove space-time correlation function and the mean- 
square displacement for a range of (small) mole fractions of 
the heavy component, and high mass ratios up to infinity. The 
latter case was realized by fixing the heavy particles at their 
starting positions during the whole simulation run. 
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